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A half-filled-band Hubbard model on an anisotropic triangular lattice {t in two bond direc- 
tions and t' in the other) is studied using an optimization variational Monte Carlo method, to 
consider a Mott transition and superconductivity arising in k-(BEDT-TTF)2X. Adopting wave 
functions with doublon-holon binding factors, we reveal that a first-order Mott (conductor- 
to-nonmagnetic insulator) transition takes place at U = f/c, approximately of the bandwidth, 
for a wide t' /t range. This transition is not directly connected to magnetism. Robust d-wave 
superconductivity appears in a restricted parameter range: immediately below Uc and the mod- 
erate strength of frustration (0.4 ^ t' /t ^ 0.7), where short-range antiferromagnetic correlation 
sufficiently develops but does not come to a long-range order. The relevance to experiments is 
also discussed. 

KEYWORDS: Hubbard model, anisotropic triangular lattice, variational Monte Carlo method, supercon- 
ductivity, d^2_y2 wave, antiferromagnetic correlation, Mott transition, fi;-(BEDT-TTF)2X 
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1. Introduction 

Organic layered compounds made of BEDT-TTF (ET) 
molecules present various characteristics of correlated 
electron systems. Among them, K-type compounds, ^"'^ 
k;-(ET)2X (X denotes an anion), have intriguing proper- 
ties as follows: (1) They have good two-dimensionality in 
conductivity with a frustrated lattice structure. (2) Most 
of them show unconventional (probably d2,2_j,2-wave) su- 
perconductivity (SC) as high as 10 K, with pseudo-gap 
behavior for T > T^. These two points are closely akin 
to high-Tc cuprates. (3) A series of k-(ET) salts give 
rise to superconductor-to-antiferromagnetic (AF) insu- 
lator transitions, through the chemical substitution of 
X, or applied pressure. (4) It was recently found that a 
compound [X=Cu2(CN)3] preserves an insulating state 
without magnetic order down to a quite low temperature 
(23 mK).'^ 

About a decade ago. Kino and Fukuyama^ studied 
the electronic structure of ET salts, applying a Hartree- 
Fock approximation to four-band Hubbard models, and 
found that if the dimerization between ET molecules is 
sufficiently strong, the conducting plane of K-type com- 
pounds is described well by a single-band Hubbard model 
on an anisotropic triangular lattice (the hopping integral 
is t in two bond directions and t' in the other) with the 
electron density of half filling. The Hartree-Fock approx- 
imation can describe a metal-to- AF insulator transition. 
Afterward, this model has been studied with various the- 
ories. Fluctuation exchange approximation (FLEX)^'^ 
and a quantum Monte Carlo method^ have yielded a re- 
sult that SC with a 1^2:2 _j^2 -wave symmetry is realized due 
to AF spin correlation (or fluctuation), as in the case of 
cuprates. Actually, experimental studies using NMR,^"^^ 
thermal conductivity,^^ and tunneling spectroscopy mea- 
surements^'^"^^ have found unconventional SC with nodes 



of the gap in k-(ET)2X. 

However, it is difficult to deal with essential aspects 
of a strong correlation, particularly Mott transitions"'^^ 
by Hartree-Fock or FLEX approximation. In contrast, 
dynamical cluster approximation (DCA)^^ and path in- 
tegral renornialization group calculation (PIRG)^* have 
shown first-order metal-to-nonmagnetic-insulator tran- 
sitions for intermediately and highly frustrated cases. 
However, both DCA and PIRG have not reached a defi- 
nite conclusion regarding SC in the vicinity of Mott tran- 
sitions. Thus, it is important to clarify the interrelation- 
ship between SC and a Mott transition. 

The purpose of this paper is to understand simulta- 
neously SC and a Mott transition arising in the Hub- 
bard model on an anisotropic triangular lattice at half 
filling, with K-ET salts in mind. To this end, we use a 
variational Monte Carlo (VMC) method, which is very 
useful in studying the physics of intermediate correlation 
strength as a continuous function of U/t}^ In the vari- 
ational wave functions, we introduce a factor controlling 
the binding between a doublon and a holon, which is nec- 
essary to appropriately treat the Hubbard model includ- 
ing Mott transitions.^" We carefully check the system- 
size dependence with highly accurate VMC data; this 
check is indispensable for the study of phase transitions. 

It is found that a first-order Mott transition takes place 
at U = Uc approximately of the bandwidth for a wide 
t' /t range. The d-wave superconducting (SC) correlation 
function indicates that robust SC is realized immediately 
below Uc for a moderate frustration (0.4 ^ t' /t ^ 0.7). 
For weakly frustrated cases {t' /t ^ 0.4), the SC is de- 
feated by an AF state with a long-range order; for highly 
frustrated cases {t' /t ^ 0.8), the SC correlation becomes 
weak for arbitrary values of U/t, together with a decrease 
in the degree of AF correlation. Thus, the origin of the 
SC can be attributed to the AF correlation. 
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Fig. 1. Lattice structure and hopping integrals of models (a) used 
in this study for k-ET salts, and (b) often used for high-Tc 
cuprates. The latter is shown for comparison. 
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Fig. 2. (Color online) Schematic explanation of a Mott transition 
due to binding of a doublon to a holon described by '^q- The 
solid (open) circles denote doublons (holons). (a) In a conductive 
regime, a doublon and a holon (charge carriers) can move inde- 
pendently, (b) In an insulating regime, a doublon and a holon 
are confined within adjacent sites; thereby, charge cannot flow. 
The results of Mott transitions are treated in S4. 



The organization of this paper is as foUows: In §2, the 
model and method used are introduced. Sections 3-5 are 
assigned to the description of our results and compar- 
isons with other theories. Section 3 treats the behavior 
of condensation energies of the d-wave and AF states. In 
§4, we demonstrate that a Mott transition actually takes 
place in the present approach by studying various quan- 
tities, and reveal the properties of this Mott transition. 
In §5, we describe robust c?-wave SC and its relevance 
to the AF spin correlation. §6 presents further discus- 
sions; first we compare our results with those of related 
experiments on k-ET salts, followed by an analysis of 
the Mott insulating state in light of strong-coupling the- 
ories. In §7, we summarize our results. In Appendix A, 
we point out the limitation of Gutzwiller-type approxi- 
mations in dealing with a Mott transition. In Appendix 
B, we comment on the related previous variational stud- 
ies in the following context. Although a couple of studies 
in the early 90's^^'^^ concluded that a Mott transition 
never arises in such wave functions as treated in this pa- 
per, this conclusion has already been proven incorrect. 
However, a recent VMC study on the present subject^'^ 
has still followed the above wrong conclusion and made 
some misinterpretations. 

Part of this study has been reported in a proceed- 

94 

mgs. 

2. Model and Method 

According to the extended Hiickel calculations based 
on the structure of k-ET salts, '^■^^ their complicated con- 
ducting planes are described with an anisotropic trian- 
gular lattice, whose unit cell is composed of a dimer of 
ET molecules. In each dimer (or unit cell), the strong 
hybridization between two ET molecules leads to a large 
split of the levels of bonding and antibonding orbitals. 
The antibonding orbital accommodates one electron, and 
composes a conduction band of half filling. Because the 
Coulomb interaction is relatively strong due to the slight 
overlap between tt orbitals, a Hubbard model has been 
considered as a realistic model. ^'^^ 

We hence consider a single-band Hubbard model on 
an anisotropic triangular lattice [see Fig. 1(a)] given as 

kcr i 

= — 2t(cos kx + cos ky) — 2t' cos{kx + ky ), (2) 



with U, t, t' > 0.^^ Throughout this paper, we fix the elec- 
tron density at half filling (n — N^/Ns — 1; A^o: electron 
number, Ns'. site number). 

To this model, we apply an optimization VMC 
method^* that can correctly treat the local correlation 
in the whole parameter space spanned by U/t and t' /t. 
As a variational wave function, we use a Jastrow type: 
'5 = V^, in which $ is a one-body (Hartree-Fock) part 
expressed as a Slater determinant, and P is a many-body 
correlation factor. Concerning V, the onsite (Gutzwiller) 
factor^^ 

= n [1 - (1 - 5)"»T«U] : (3) 
i 

was found effective for i- J-type models. For the Hub- 
bard model, however, it has been known'^^ that the wave 
functions Vg^{= '^g) yield several unfavorable results, 
especially for large values of U/t: (1) The momentum dis- 
tribution function n(k) is an increasing function of |k|. 
(2) 2k-p anomalies in the spin [charge density] structure 
factor S'(q) [-/V(q)] cannot be properly represented. (3) 
A Mott transition cannot be described, in addition to a 
sizably high variational energy. 

To overcome these flaws, one has to introduce an inter- 
site correlation factor. ^^■^'^ In particular, at half filling, 
the doublon-holon binding factor Vq,"^^'^^'^* in addition 
to Vg , is an indispensable element for describing a Mott 
transition. Generally, Vq should be written as a func- 
tional of /i(r), which is a function of the relative position 
of a doubly occupied site (doublon) from an empty site 
(holon), r.^^ Here, however, we integrate only the two 
major elements into Vq, namely, for nearest neighbors, 
fi [r = {±x, 0) and (0, ±y)] and for second (diagonal) 
neighbors, ^' [r = {x,y) and (— x,— j/)]. 

^Q = n(i-^QO(i-M'Qr') (4) 

i 

Ql^^'^ = n [^^^(1 - e,+rir')) + e.(l - d,+rir'))] (5) 
r(r') 

Here, di — rii-^nii, Ci = (1 — ni|)(l — riij), and t (r') runs 
over all adjacent sites in the bond directions oft (t'). Our 
trial wave function thus becomes 

*Q - VQVG<i>. (6) 

The factor Vq is naturally derived by taking account 
of virtual states in the second-order processes in the 
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strong- coupling expansion, ^^'^^ and controls the binding 
strength between a doublon and a holon with the param- 
eters ji and ^' (/z, /z' < 1). For ~ 0, doublons and 
holons, which are charged, move around almost freely, 
whereas for 1, a doublon and a holon are tightly 
bound within adjacent sites in the bond directions of t 
{t') (See Fig. 2). Thus, current cannot flow for /i — > 1. 
It has been argued in detaiP° that this factor Pq can 
describe a Mott transition and its existence has been ac- 
tually confirmed for a regular square lattice. Afterward, 
Mott transitions have been found using Pq for various 
systems (refer to §B.l in Appendix B).^^'^'^"^'^ 

Regarding <I>, wc treat three cases in this paper, a d- 
wave BCS state, an AF state and a Fermi sea. First, we 
introduce a singlet wave function of the form of the BCS 
function with a dx2-y2-waye gap 



l^¥'k4^clkj |0), 



ek-C + \/(£k-C)2+A^' 



(7) 



(8) 



with 



Ek = — 2t(cos kx + cos ky) — 2i' cos{kx + ky), (9) 

and Ak = A£i(cos fcj; — cos ky). Here, (, Ad, and i' are var- 
iational parameters to be optimized and ( corresponds to 
the chemical potential. Although the parameter indi- 
cates the magnitude of the d-wave gap, it does not nec- 
essarily mean the stability of SC; in the insulating phase, 
Ad is related to the spin gap. Furthermore, we take ac- 
count of a renormalization effect of the quasi-Fermi sur- 
face due to the electron correlation. The main part of 
this effect is introduced into $d by optimizing i' /t, in- 
dependently of t'/t given in the Hamiltonian eq. (1).'^^ 
As we will sec later, t'/t notably deviates from the orig- 
inal t'/t, particularly in the insulating phase, and the 
renormalization affects the stability of the d-wave singlet 
state. 

Next, as a competing state for weakly frustrated 
(small-t'/t) cases, we consider the commensurate AF 
long-range ordered state 



^AF = n (^kcL + Sgn((7)i}k4+K<T) |0) ' 



k,cr 



with 



\ 



(+) 



£k 



v4 + Aip^ 



(10) 



(11) 



Here, K = (tt, tt), and the wave- number space (k) is re- 
stricted to the folded AF Brillouin zone. Aaf is a var- 
iational parameter indicating the magnitude of the AF 
gap. At half fiUing, *af(= PqPg^af) with finite Aaf 
is always insulating, as in the mean-field theory. 

Lastly, as a reference state, we use a projected Fermi 



t 

k| </cfO" '""'ker 



Al- 



sea: ^f{= PqPg^f) with *f = O 
though \1/f also brings about a Mott transition at a some- 
what larger U than the bandwidth, we will not treat 



the insulating regime of ^'p in this paper, and thus call 
it the 'normal' or 'metallic' state. 

In the remainder of this section, we briefly mention the 
details of VMC calculations. Because the trial functions 
we treat have at most six parameters, we employ a sim- 
ple linear optimization of each parameter with the other 
parameters fixed. In a round of iteration, every param- 
eter is once optimized in one dimension. In almost all 
the cases of optimization in this study, the parameters 
converge within five rounds, after which we continue the 
optimization for more than ten rounds. The optimized 
values of the parameters and energy are determined by 
averaging the results of the ten rounds after convergence. 
Because the optimization is carried out with 2-5 x 10^ 
samples, our data are substantially averages of several 
million samples. Thus, we can greatly suppress the sta- 
tistical errors in energy down to the order of 10~H. Such 
precision is necessary to correctly analyze the physics in 
the vicinity of a phase transition. Physical quantities are 
calculated with the optimized parameters thus obtained 
with 2.5 X 10^ samples. 

The systems used are of Ng = LxL sites with periodic- 
antiperiodic boundary conditions. A closed shell condi- 
tion is always imposed. We mainly study the systems for 
L — 10 and 12, and check those for L = 14 and 16 when 
we would like to know the system-size dependence more 
definitely. As we will see later, for t' ^ 0, in contrast to 
for t' = 0, the system-size dependence is not monotonic; 
in particular, the SC correlation function is very sensitive 
to the configuration of occupied k points. Thus, we need 
to deduce the properties for L ^ oo from the tendency 
of finite-sized systems in varying L. 

3. Stability of d-wave and AF states 

In this and the following two sections, we describe re- 
sults of VMC calculations, and compare with other theo- 
retical ones. In this section, we first consider the compe- 
tition between d-wave and AF states, and then mention 
the critical behavior found for '^q. 

To consider the competition between and we 
take up the condensation energies defined by E^{E^^) = 
- E'^{E^'^), in which E'^, E^'^ , and E^ are the op- 
timized variational energies per site with respect to ^q, 
\E'q^ and , respectively. For t' = (the regular square 
lattice), it is well-known that the ground state has an 
AF long-range order for arbitrary positive values of U/t 
due to the complete nesting condition; this feature can 
be properly described by the present approach, namely 
E^^ > E^ > 0.19 The maximum E^^ is O.USt {U/t = 8) 
for L = 12, and will be a little larger for L oo."^^ 

In Fig. 3, we show E^ and E^^ for four values of t'/t 
(= 0.2-0.8). As the frustration becomes stronger {t'/t 
increases), the magnitude of E^^ rapidly decreases, and 
thereby the stable region of the AF state shrinks. For 
t'/t = 0.2, the AF region still remains: U ^ 6.7 {L = 10) 
and U ^ 7.5 {L = 12). For t'/t = 0.4, however, the AF 
region becomes very narrow for L = 10 (5 ^ t/ ^ 6.5), 
and finally vanishes for L = 12 {E^^ < E^), as seen in 
the inset of Fig. 3(b). These indicate that at t'/t ^ 0.4, 
the ground state switches from the AF state to the d- 
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Fig. 3. (Color online) Condensation energies E^jt of d-wave and AF states for four values of t' jt: (a) t' jt = 0.2, (b) 0.4, (c) 0.6, and 
(d) 0.8. For t' /t = 0.6 and 0.8, is zero for arbitrary values of U/t. The arrow in each panel indicates the critical point, Uc/t of 

the Mott transition arising in the d-wave state. The values of Uc/t shown are estimated from the quantities treated in §4 and §5. The 
inset in (b) is a close-up near the critical point. Systems of mainly L = 10-12 are compared. 



wave state, although the system-size dependence of E^^ 
is not simple near t'/t = 0.4 (see rus in Fig. 13). For t'/t = 
0.6 and 0.8, the AF state is never stabilized with respect 
to the normal state. This feature is natural at half filling, 
because the nesting of the Fermi surface monotonically 
becomes worse, as t'/t increases. 

In the following, we consider the behavior of E'^. For 
small values of U/t { ^ A), E^ is very small for arbi- 
trary values of t'/t, indicating that SC is fragile, if any. 
E^ gradually increases for intermediate values of U/t: 
A ^ U/t <, 6.5 for t'/t = 0.2, 5.5 <, U/t <, 6.5 for 
t'/t = 0.4, and 6 ^ C//i ^ 7 for t'/t = 0.6. For t'/t = 0.8, 
this gradual increase is hardly seen. As will be shown in 
§5, this increase corresponds to the stabilization due to 
robust SC. For every t'/t, E^ suddenly and markedly in- 
creases Situ — Uc, which is indicated by arrows in Fig. 3. 
This increase in E'^ certainly indicates the stabilization 
of the d-wave singlet state, but is not necessarily a 
SC state, as explained in §2. We will pursue this behavior 
in the next section. 

4. Mott transitions 

In this section, we first demonstrate that the critical 
behavior found in in §3 is attributed to a Mott tran- 
sition by studying various quantities (§4.1), and then re- 
veal the properties of the transition from different points 
of view (§4.2). 



^.1 Confirmation of Mott transitions 

First, we consider the optimized values of the vari- 
ational parameters in \E'q, which are plotted in Fig. 4 
for four values of t'/t. For each t'/t, all the varia- 
tional parameters have discontinuities at the same U/t, 
which is identical to Uc/t mentioned in the preced- 
ing section. Here, Uc is determined accurately: Uc/t — 
6.55,6.75,7.15, and 8.05 for t'/t = 0.2, 0.4, 0.6, and 
0.8, respectively. These discontinuities indicate that some 
first order phase transition takes place At U = Uc- Actu- 
ally, near U = Uc, E'^ has double minima in the parame- 
ter space; a hysteresis in E''' is observed, as to whether we 
approach Uc from the small-J7 side or the large-J7 side, 
and particularly notable in the cases of large t'/t. Here, 
recall that the doublon-holon binding factor /i is a good 
index of a Mott transition. In the small-C/ regime, /i is 
as small as 0.1-0.3 and almost independent of L, whereas 

is likely to approach 1 ior U > Uc, in allowing for the 
large L dependence [Fig. 4(d)]. Thus, this transition is 
expected to be a Mott transition. In the following, we 
will confirm this expectation with various quantities. 

Let us start with the momentum distribution function 

"(k) = ^E(4.ck.). (12) 

Shown in Fig. 5 is n{\<L) measured with the optimized 
parameters along the path (0, 0)-(7r, 0)-(7r, 7r)-(0, 0) for 
t'/t = 0.4 and 0.8. Because vPg has a node of the gap 
in the (0, 0)-(7r, tt) direction, metallic properties should 
develop in this direction; n(k) should have a discontinu- 
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Fig. 4. (Color online) Optimized values of variational parameters in d-wavc singlet state, ^g, for four values of t' /t as functions of U/t; 
(a) g [onsite (Gutzwiller) factor], (b) t'/i [band renormalization factor], (c) A^j/i [d-wave gap], (d) /i [doublon-holon binding factor in 
the direction of t], (e) ^' [the same of t'], (e) C, [chemical potential]. The systems are of L = 10-14. Statistical errors in i' /t and in C, 
increase as U/t decreases or t' /t increases, due to the discrete energy levels or k points. 



ity at fcp. For both values of i'/i, the discontinuity at fcp 
is apparent for U < Uc, whereas, for U > Uc, the discon- 
tinuity, that is, the quasi-Fermi surface disappears. Thus, 
metalhc properties are abruptly lost at U — Uc even in 
the nodal direction of the dj.2_y2 wave. 

To see this feature quantitatively, the quasi-particle 
renormalization factor Z is a suitable indicator, which 
roughly corresponds to the inverse of effective mass, un- 
less the /c-dependent renormalization of self energy is se- 
vere. In Fig. 6, we plot Z estimated from the magnitude 
of jump in n(k) at k = kp in the nodal direction."*^ For 
all the values of t'/t, Z decreases slowly for U < Uc, 
whereas, at U ^ Uc, Z suddenly vanishes with a siz- 
able discontinuity, reflecting the first-order character of 
the transition. The system-size dependence of Z is weak, 
and the magnitude of Z tends to increase as L increases 
{t'/t — 0). The behavior of Z strongly suggests that the 
effective electron mass diverges for U > Uc- Incidentally, 
the discontinuity in Z at C/ = J7c is a feature that dif- 



fers from those of ^'g^^' and of a dynamical mean-field 
theory for the hypercubic lattice, '^'^ in which Z gradu- 
ally decreases and vanishes aX U = Uc without a jump. 
Thus, the first-order character of the transition is more 
conspicuous in ^g. 

Next, we consider the charge structure factor 

^(q) = e^^<^^~^^^ {N,N,) - n\ (13) 

with Ni — ni-^ + Within the variational theory, 
A^(q) cx |q| for |q| — > 0, if the state does not have a gap 
in the charge degree of freedom, whereas A^(q) oc q^, 
if a charge gap opens. In Fig. 7, A^(q) is depicted for 
two different values of t'/t (0.4 and 0.8) near U = Uc- 
For both values of t'/t, N{q) near the F point (0,0) 
seems linear in |q| for U < Uc (solid symbols), whereas 
it abruptly changes to being roughly quadratic in |q| for 
U > Uc (open symbols). It follows that \E'q is gapless 
in the charge sector and is conductive for U < Uc, but a 
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Fig. 5. Momentum distribution function n(k) of d-wave SC state 
for (a) t'/t = 0.4 {Uc/t ~ 6.75) and (b) 0.8 {Uc/t ~ 8.05). The 
sohd (open) symbols denote the points oi U < Uc {U > Uc). 
The arrows indicate the positions of quasi-Fermi surface in the 
node-of-gap direction (0,0)-(7r, tt). The systems are of L = 12. 
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Fig. 7. Charge structure factor N{q) of d-wave singlet state for 
(a) t'/t = 0.4 and (b) 0.8. Solid (open) symbols indicate the data 
for U <Uc {U > Uc). The systems are of L = 12. 
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Fig. 6. (Color online) Quasiparticle renormalization factor Z in 
node-of-gap direction of d-wave singlet states for four values of 
t'/i as a function of U/t. The deviation of the data for t' /t = 0.4, 
L = 12 and U/t = 6.5 is a special case, and is ascribed to the 
discrete k points of the system."*^ 



charge gap opens for U > Uc and ^'g becomes insulating. 

As further evidence for a Mott transition, we consider 
the doublon density (density of doubly-occupied sites), 

1 (Hint) 



Ns U 



(14) 



which is regarded as the order parameter of Mott transi- 
tions,^'* by analogy with the particle density in gas-liquid 
transitions. In Fig. 8, we plot D for several values of t'/t. 
As the interaction strength U/t increases, the magnitude 
of D decreases almost linearly in U/t for small U/t, but 
abruptly drops at Uc/t, and then decreases slowly for 



0.15fe 




0.05 



U/t 

Fig. 8. (Color online) Density of doublon (doubly occupied site) 
of d-wave singlet state as a function of U /t near U = Uc for four 
values of t'/t. For t'/t = 0, the horizontal axis is shifted by 0.5 
(see upper axis) for clarity. This quantity is reduced to 0.25 (0) 
in the limit of U/t — » (oo) at half filling. Data for some system 
sizes are shown for each t'/t. 



U > C/c-'*^ This behavior of D corroborates the idea of a 
first-order Mott transition. 

Owing to the behavior of all the quantities studied 
above, we can safely judge that a Mott transition takes 
place in the d-wave singlet state at U = Uc for arbitrary 
values of t'/t. 

4-2 Properties of Mott transitions 

For a start, we contrast this transition with the well- 
known theory of Brinkman and Rice,**^ which is based 
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on '^Q. Their theory is a milestone in the research on 
Mott transitions, and stiU survives as a symbohc concept 
of a certain class of transitions. However, it should be 
noted that the result of their approximation additionally 
applied, the so-called Gutzwiller approximation^^ (GA), 
has been proven correct only in infinite space dimen- 
sions. In accurate calculations of \I>Qj20,32,47,48 ^ nietal- 
insulator transition is absent in finite dimensions. We 
sum up the properties of the Brinkman-Rice transition. 
For a half-filled-band Hubbard model, as U/t increases 
toward U = Ubr = 8Eq {Eq: energy for {7 = 0), elec- 
trons are gradually localized in lattice sites one by one, 
and the magnetic susceptibility diverges as C/ — s- Ubr', 
this transition is a continuous type. For U > Ubr, elec- 
trons are completely localized and never move, namely 
E^n = Eint = 0. Hence, the behavior of \E\ ^ At^ /U 
(= J) in the insulating regime cannot be derived from 
this framework. 

In contrast, in the present theory, the Mott transition 
is described not by the disappearance of charge carriers, 
but by the binding (and unbinding) of negative carri- 
ers (doublons) to positive carriers (holons), as explained 
for Vq in §2 and schematically shown in Fig. 2. In this 
mechanism, doublons (and holons) never vanish even for 
U > ?7c, as we can actually observe 5 > in Fig. 4(a) and 
D > in Fig. 8. The process of creating a temporary pair 
of a doublon and a holon through the so-called virtual 
hopping leads to the behavior proportional to t^/J7. 

In this connection, a couple of studies^^'^^ have re- 
cently addressed the present subject, applying GA for 
d-wave BCS states^^ to Hubbard-type models with ex- 
tra exchange-coupling terms. We compare the results of 
these studies with those of conventional GA for the Fermi 
sea in Appendix A. As long as the Mott transition is con- 
cerned, the results of GA for the c?-wave BCS state have 
features qualitatively different from the present VMC re- 
sults, and share the behavior with the conventional GA 
for the Fermi sea. 

Next, we consider the function of ^' . For U < Uc, 
/i' gradually increases, as U/t increases, similarly to fi, 
although l^'l is small compared with \^\ [Figs. 4(d) and 
4(e)]. Thus, /Li' plays an analogous role to /i in the conduc- 
tive phase. Concerning the t' /t dependence, the behavior 
of fi' is opposite to that of /i, as naturally expected. At 
U — Uc, however, /i' suddenly drops to a negative value, 
irrespective of t' /t. This behavior means that, in the di- 
agonal direction, a doublon and a holon as two indepen- 
dent particles are more favorable than those as a bound 
pair in the neighbors. As we will demonstrate shortly, 
this is probably because the AF short-range correlation 
noticeably develops for U ^ Uc', thereby, parallel spin 
configurations are favored in the diagonal directions, and 
suppress the hopping in terms of t' . This marked direc- 
tion dependence of /i(r) is a new feature specific to two 
dimensions, in particular, to the d-wave singlet state. Ac- 
cording to the exact-diagonalization calculations for one- 
dimcnsional systems, the behavior of ^(r) for a large U /t 
is simple, namely, ^(r) is approximately proportional to 
exp(— r/ro), in which ro is a constant. 

In what follows, we consider the t' /t dependence of 
this Mott transition in some points. In Fig. 9, the crit- 



9 





W/t 










> 




rf-wave 



^0 0.2 0.4 0.6 0.8 1 
t'/t 



Fig. 9. (Color online) Comparison between bare bandwidth W 
and critical value Uc of Mott transition for d-wave singlet state. 
Because Uc/t is determined using finite-sized data (mainly L = 
12), the critical values for L = 00 may somewhat shift to the 
larger side of U/t, especially for large values of t'/t. 

ical value Uc is compared with the bare bandwidth W. 
For the present dispersion eq. (2), W is constant, 8t, for 
< t'/t < 0.5, and gradually increases as t'/t increases 
over 0.5. As W/t increases, the absolute value of kinetic 
energy increases, leading to a situation advantageous for 
the conductive state. Consequently, Uc/t increases with 
W/t; actually, Uc/t starts to increase at a smaller t'/t 
than W/t. This feature is common to different Mott- 
transition systems. The ratio Uc/W has become a topic 
in some studies; both dynamical mean field calculations 
for infinite-dimensional lattices, ^'^ and VMC calculations 
using \l/^20,37 Jqj. ^/^^ _ Q yigi(j values of Uc/W some- 
what larger than 1. In this work, however, Uc/W = 0.81 
for t'/t — 0, and increases with t'/t toward 1.02 (L — 12) 
for t'/t = 1.0, which is relatively close to the result of 
DCA, 1.15-1.2, for t'/t = 1." Conversely, this ratio ob- 
tained using PIRG is rather small, Uc/W — 0.58, even 
for t'/t = 1.0.18 

Now, we turn to the order of the Mott transition. Re- 
turning to Fig. 4, we notice that the discontinuity at 
U — Uc increases, as t'/t increases, in every parameter. 
Such a tendency is common to other quantities like Z 
(Fig. 6) and D (Fig. 8), and is also observed in the hys- 
teresis or double-minimum structure in E^, which be- 
comes more conspicuous, as t'/t increases (not shown). 
This indicates that the character of a first-order tran- 
sition becomes prominent, as the frustration becomes 
stronger (t'/t increases). This feature seems unique to 
and in contrast to that of the Mott transition in 
as well as that of PIRG.'^^ In these cases, the ef- 
fect of frustration makes the transitions more like those 
of the continuous type. This characteristic behavior of 
is probably caused by the prominent stability in the 
insulating regime, studied in the following and §6.2. 

In the remainder of this subsection, we consider the 
spin degree of freedom of the Mott insulating state de- 
scribed by ^'q for U > Uc- Shown in Fig. 10 is the spin 
structure factor, 

Si^)-^^J2^'''''-'''''''HS!S]), (15) 

for two values of t'/t. Leaving the behavior in the con- 
ductive regime for the next section, here we focus on the 
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Fig. 10. (Color online) Spin structure factor 5(q) of d-wave sin- 
glet state along the same path as in Fig. 5 are compared between 
(a) t' /t = 0.4 and (b) 0.8. Solid (open) symbols denote the points 
belonging to the conductive (insulating) region. The insets in 
each panel arc close-ups near the F point (0, 0) in the directions 
of (0, 0)-(7r, 0) (left) and (0, 0)-(7r, vr) (right). The systems are of 
L = 12. 



insulating state shown by open symbols. For U > Uc, 
S'(q) has a prominent peak at q = K = (tt, tt), compared 
with that for U <Uc\ the AF spin correlation markedly 
develops in the insulating phase. On the other hand, the 
height of peaks, S'(K), has only weak system-size depen- 
dence compared with L as well as L^: For example for 
t'/t = 0.4 and U/t ^ 7, S(K) = 6.36, 6.85, and 7.28 for 
L — 10, 12, and 14, respectively. In the same sense, the 
staggered magnetization 



1 



(16) 



with = (riii — nii)/2, is always zero for "ifg within 
the statistical fluctuation. Thus, this insulating state has 
sizable short-range AF correlation, although a long-range 
order is not formed. The enhancement of AF correlation 
reflects the change in some quantities. In Fig. 4(b), where 
the renormalized value of the diagonal hopping integral 
is plotted versus U/t, we find that t'/t suddenly drops to 
very small values at [/ = Uc, irrespective of the model 
value t'/t. Therefore, the quasi- Fermi surface recovers the 
nesting condition to a considerable extent, as shown in 
Fig. 11(c). This feature was also reported in a study us- 
ing DCA.^^ Moreover, the chemical potential ^ becomes 
almost zero [Fig. 4(f)], namely, the value of the regular 
square lattice (<' — 0), for which the nesting condition is 
fully satisfied. 

The short-range nature of the AF spin correlation 
affects the low-energy spin excitation. In the insets of 
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Fig. 11 



Effective quasi-Fermi surfaces for (a) U/t = 0, (b) 6 (SC 
phase) and (c) 8 (insulating phase). Here, the model value t'/t 
is fixed at 0.4, whereas the optimized value of t'/t in the wave 
function is renormalized as i' /t = 0.31 (0.02) for U/t = 6 (8). 



Fig. 10, we depict the magnifications of S{q) near the F 
point. For U > Uc, S{q) is not proportional to |q| for 
|q| but to |qp in both bond and diagonal directions 
and in both values of t'/t, indicating that this insulating 
state has a gap in the spin degree of freedom. 

In §6.2, this insulating state will be considered again 
in light of strong coupling theories. 

5. Superconductivity and phase diagram 

In this section, we elucidate where, in the parameter 
space, the c?-wave SC is dominant (§5.1). Therefore, we 
construct a phase diagram of the present model. We also 
study the properties of the SC, emphasizing the relevance 
to the AF spin correlation (§5.2). 

5.1 Range of robust superconductivity 

As an indicator of SC with the d^^-y^ symmetry, it is 
convenient for the present approach to study the d-wave 
nearest-neighbor pair correlation function^^ defined as 



1 

4iv: 



i t,t'— x.y 



(-1) 



1-(S(t,t') 



(At(R,)A,,(R,-f r)). 



(17) 



in which x and y denote the lattice vectors in the x- 
and y-directions, respectively, and A!jl(Ri) is the creation 
operator of a nearest-neighbor singlet expressed as 



(18) 



If Pd(r) has a finite value for |r| — > oo, one can judge 
that an off-diagonal long-range order is realized. How- 
ever, for finite systems, especially, in the cases with small 
U/t, where the correlation length is large, we have to 
measure long-distance values of ^^(r) appropriately. In 
this paper, we average Pd(r)'s with r being on the line 
segments (0, L/2)-(L, L/2) and (L/2, 0)-(L/2, L), which 
correspond to the furthermost edges from the origin, 
for weak correlation regimes.^'' For strong correlation 
regimes, i-'rf(r) becomes almost constant for |r| > 3;^^ 
we average Pd(r)'s with |r| > 3 in these cases. 

In Fig. 12, we plot the averages of Prf(r) thus obtained 
(P|™) as a function of U/t for four values of t'/t. The 
data for t'/t = 0.7 are plotted in Fig. B-2(b). First, let 
us consider features common to the cases of weak and 
intermediate frustrations {t'/t = 0.2-0.7). For small U/t, 
the magnitude of is as small as that for the non- 
interacting case {U/t ~ 0); here, robust SC cannot be 
expected. P^™ starts to increase appreciably at a value 
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Fig. 12. (Color online) Real-space SC correlation function for nearest-neighbor pairs of d-wave symmetry for (a) t'/t = 0.2, (b) 0.4, 
(c) 0.6, and (d) 0.8. Plotted are the averages of the long-distance part of PdiJ^) explained in the text. The insets in (b) and (d) are 
close-ups near Uc/t. Data of different system sizes (L = 10-14) are compared. The same quantity for t' /t = 0.7 is shown in Fig. B-2(b) 
in Appendix. 



of U/t (Uonsct/t), which is different for a different t'/t. 
There is a peak in PJ'™ immediately below Uc- The fact 
that the magnitude of PJ'™ near the peak is almost inde- 
pendent of system size indicates that a region of robust 
SC exists immediately below the Mott critical point Uc- 
At Uc/t, P|™ suddenly drops to a very small value, which 
is justifiably expected from the Mott transition. In the 
insulating regime {U > Uc), as the system size L in- 
creases, the magnitude of -P^™ rapidly decreases as seen 
in every panel in Fig. 12, and probably vanishes in the 
limit of L — > oo. In the insulating regime, the statistical 
fluctuation in the VMC data is much smaller than that 
in the conductive regime; the disappearance of Pd{r) for 
U > Uc is absolutely certain. In §B.2, we comment on 
a recent VMC study,^'^ which has drawn contradictory 
conclusions from these results. 

Next, we proceed to the t'/t-dependence of P|™. 
As the strength of frustration t'/t increases, the onset 
Uonsct/t, at which P_^™ starts to increase appreciably, 
increases: Uonset/t 4 (t'/t = 0.2), - 5 (0.4), ~ 6 
(0.6), and ~ 7 (0.7). Because, compared with the in- 
creasing rate of Uonsct/t, that of Uc/t is slow: 6.55 (0.2), 
6.75 (0.4), 7.15 (0.6), and 7.45 (0.7), the range in U/t 
of the dominant SC becomes narrower, as t'/t increases, 
as observed in Fig. 12. In Fig. 13, we depict the max- 
imal PJ'™ with t'/t fixed and U/t varied. This value 
is almost constant for t'/t < 0.4, decreases as t'/t in- 
creases beyond t'/t = 0.4, and finally vanishes at ap- 
proximately t'/t = 0.8. Actually for a strongly frustrated 
case {t'/t = 0.8), an enhancement in Pj'™ is hardly seen 




0.4 0.6 
t'/t 

Fig. 13. (Color online) Comparison between the maximal values 
of P|^° (average of the SC correlation function) for fixed values 
of t'/t (with U/t changed) and S(K) (spin structure factor at the 
AF wave number) calculated with the same parameter of P^^" ■ 
Also shown is the sublattice magnetization rris calculated using 
at U/t = 6.0. The scatter of data points, in particular, 
near the critical values of t'/t, is ascribed to finite system sizes. 
Statistical errors are much smaller. 



even below Uc/t (=8.05), as in Fig. 12(d), in which the 
magnitude of PJ'™ obtained for U/t — 4-8 is compa- 
rable to the statistical fluctuation. Thus, stable SC is 
absent from the strongly frustrated (or nearly isotropic) 
region: t'/t ^ 0.8. Now, recall that the AF state with 
a long-range order is more stable than the SC state for 
t'/t ^ 0.4, as demonstrated in §3. This feature is reflected 
in the sublattice magnetization mg [eq. (16)] for , as 
shown in Fig. 13. Thus, we estimate the range of stable 
SC as 0.4 <, t'/t <, 0.7. 
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Incidentally, as shown in the insets of Fig. 10(a), 5(q) 
tends to be proportional to |qp for |q| — > as U /t in- 
creases ([/ < Uc), indicating a spin gap is apt to open as 
SC develops. In contrast, in the insets of Fig. 10(b), S'(q) 
remains linear in |q|, even if U/t increases; remains 
metallic for t' /t — 0.8, even if U approaches Uc- 

The present results for SC are consistent, even quan- 
titatively, with those of FLEX.^'^ According to Kino 
and Kontani, the AF order defeats SC for t' /t < 0.4, 
whereas the region of finite of (i^2_j^2-wave SC appears 
for 0.4 < t'/t < 0.8 with the highest values taken for 
t'/t = 0.5-0.6. Finally, of SC disappears for t'/t > 0.8. 
Furthermore, the range of SC in U/t displays a tendency 
similar to that of the present study. 

Having discussed the range of SC in the t'-U space and 
its stability, let us now reconsider how various quantities 
behave in this SC region. First, concerning the conden- 
sation energy shown in Fig. 3, we are now aware that 
the gradual increase in E^. for J/onsot ^ U < Uc derives 
from SC, whereas the sizable gain for U > Uc is not 
due to SC but to the d-wave singlet spin state. Simi- 
lar gradual changes can be seen in the variational pa- 
rameters (Fig. 4); in particular, the increase in for 
^^onsct ^ U < Uc directly corresponds to the development 
of SC [Fig. 4(c)]. The slow decrease in i' /t [Fig. 4(b)] 
means that the degree of renormalization of £k is low in 
the SC state; the quasi-Fermi surface is not very differ- 
ent from that for U/t = as in Figs. 11(a) and (b). This 
tendency is again consistent with the result of FLEX.^ 
Note that such a gradual change never appears in any 
quantity in the case of t'/t ~ 0.8, where robust SC is 
absent. 

5. 2 Properties of superconductivity 

Because the SC develops as soon as the AF long-range 
order disappears, the AF correlation must be important 
for the mechanism of SC. We hence begin with the spin 
structure factor S'(q) [eq. (15)], returning to Fig. 10. For 
U/t = 0, S{q) does not have a peak but rather a dip at 
the AF wave number K = (vr, tt) for both t'/t — 0.4 and 
0.8, due to the frustration. In the weakly frustrated case 
{t'/t = 0.4), in which robust SC appears, S(K) steadily 
grows as U/t increases, and K becomes a characteris- 
tic wave number even in the conductive regime, U < Uc 
[solid symbols in Fig. 10(a)]. In contrast, in the strongly 
frustrated case {t'/t = 0.8), in which the SC correlation 
does not develop, although 5(q) somewhat increases as 
a whole for U < Uc [solid symbols in Fig. 10(b)], pre- 
serving the shape for C//t = 0, yet 5(q) does not exhibit 
a special enhancement at q = K, but has maxima at 
incommensurate wave numbers. To check this point, we 
plot, in Fig. 13, 5(K) for 'I'q giving the maximal PJ^™; 
when S'(K) decreases, the maximal simultaneously 
decreases. We have confirmed it for a wide range of model 
parameters that whenever -P^™ is appreciably enhanced, 
S'(q) has an evident peak at q = K. From these discus- 
sions, we can conclude that SC in this model is induced 
by the AF spin correlation. Thus, we can understand that 
the shape of the quasi-Fermi surface of '^q tends to be 
more square-lattice-like {i'/t' < 1), as mentioned earlier. 
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Fig. 14. (Color online) Phase diagram in the t'-U plane con- 
structed from the present VMC results. 

so as to place the nesting vector K on hot spots having 
a larger density of state. It is probable that the mecha- 
nism of SC proposed here is basically the same as those 
of the weak-coupling theories,^' "^'^^ and is also common 
to that of high-Tc cuprates at half filling, if any, which 
have quite similar lattice structure (see Fig. 1).'^'' 

Here, we summarize the features of this SC in light 
of the AF correlation. (1) For small values of U/t, be- 
cause magnetic correlations do not fully develop, robust 
SC does not appear. (2) For small values oft' /t, although 
the AF correlation develops as U/t increases, an AF long- 
range order defeats SC due to the nesting condition. (3) 
For moderate values of t'/t as well as U/t, the AF cor- 
relation is not strong enough to form an AF long-range 
order, but the short-range part survives. Thereby, stable 
d-wave SC is realized. (4) For large values of t'/t, even 
short-range AF correlation decays due to severe frustra- 
tion. Thus, SC is not likely to arise. (5) For considerably 
large values of U/t, a Mott insulating state suppresses 
SC. 

On the basis of the results mentioned upto this point, 
we constructed a phase diagram in the model space 
spanned by t'/t and U/t (Fig. 14). Robust d-wave SC 
appears in the highlighted area (0.4 ^ t'/t 0.7, 
C^Dnsct ^ U < Uc), which is in contact with both AF 
and Mott insulating areas, and continues to the metal- 
lic area in the small-C/ /t or large-i'/i sides, where the SC 
correlation is very weak. The conductive (metallic or SC) 
and nonmagnetic insulating areas are divided by a first- 
order Mott transition described within the d-wave singlet 
state. The boundary between the areas of long-range AF 
insulator and the nonmagnetic insulator (or supercon- 
ductor) is determined by comparing the energies of ^'g^ 
and ^Q. 

Finally, we take up another weak-coupling feature of 
this SC. Shown in Fig. 15 are the components of con- 
densation energy, that is, the differences in kinetic and 
interaction energies between and ^'g given as 

^-^kin = -^kin - ^kiny (19) 
^Ef^t = ^int-^itt> (20) 

in which AE^-^ + AEf^^^ = E^ {> 0). In the SC regime 
{U < Uc), AEf^^ (AiJ^jj^) is always positive (negative). 
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Fig. 15. (Color online) Differences in (a) interaction and (b) ki- 
netic energies between normal and d-wave states for t' jt = 0.2, 
0.4, 0.6, and 0.8. 

regardless of t' jt. This indicates that the SC transition 
is induced by the gain in interaction energy at the cost 
of kinetic energy. This feature smoothly continues to the 
weak-coupling limit {U jt — > 0), and common to conven- 
tional BCS superconductors. Although the component of 
energy gain switches to the kinetic energy for U ^ lOt, as 
seen in Fig. 15, yet a SC phase is not realized for U > [4; 
the SC due to the kinetic-energy gain is not realized here, 
in contrast to the doped cases. 

6. Further Discussions 

In §6.1, we compare the results obtained in this work 
with those of related experiments on k-ET salts. In §6.2, 
as a continuation of §4.2, we analyze for U > U^va 
light of strong-coupling theories, and propose a couple 
of improvements. 

6.1 Comparisons with experiments 

First, let us think of SC. If the present results for SC is 
applicable to k-ET salts, the SC arising in them should 
be induced by AF correlation, which becomes weak as 
t' /t increases (Fig. 13). Actually, in a class of k-ET salts, 
a compound with weaker frustration (smaller t' /t) tends 
to have higher T^: for X=CuN(CN)2Br {f /t = 0.68), 
Cu(NCS)2 (0.84), and Cu2(CN)3 (1.06), = 11.6, 10.4, 
and 3.9K (under pressure of 0.4GPa), respectively. Note, 
however, that, according to our result, robust SC does 
not arise for t' /t ^ 0.8; this fact seemingly contradicts 
the experiments on k-ET salts, particularly of the com- 
pound X=Cu2(CN)3. A point we must not overlook here 
is the effect of pressure. The compound X=Cu2(CN)3 
has an almost regular triangular structure {t' — t) under 



ambient pressure, which fact is also confirmed by a high- 
temperature-expansion study.^" However, this compound 
does not exhibit SC, unless a pressure of 0.4-0.6GPa is 
applied. It is well-known that applied pressure gen- 
erally enhances the bandwidth (reduce U/t). It is rea- 
sonable to expect that t'/t, in addition to U/t, should 
be greatly varied by pressure, if we allow for the fact 
that organic compounds are soft and anisotropic, and 
the values of t and t' are sensitive to bond angle. Ac- 
tually, a recent NMR-Ti experiment for X=Cu2(CN)36i 
has suggested the enhancement of AF fluctuation under 
pressure; in this case, it is likely that t'/t is sufficiently 
reduced to induce SC. Although the pressure effect on 
t'/t is quite an important problem, it has not yet been 
established whether t'/t is reduced^^ or enhanced^'* by 
pressure. 

Next, we touch on Mott (SC-nonmagnetic insulator) 
transition. The existence of a first-order Mott transi- 
tion is consistent with the experiment on temperature- 
dependent magnetic susceptibility and the nuclear spin- 
lattice relaxation rate of k-(ET)2Cu2(CN)3.'* In this 
compound, however, the ground state in the insulating 
phase has a low-energy spin excitation; the spin gap will 
be very small, if it exists. This point is not necessarily 
consistent with the present result. Furthermore, the in- 
sulating states of most k-ET salts exhibit a long-range 
AF order, although they are considered strongly frus- 
trated with t'/t ^ 0.6. Even if we take account of the 
fact that actual materials possess the effect of three di- 
mensions, the trial function in the Mott insulating regime 
leaves room for improvement with regard to the spin cor- 
relation, as will be mentioned in the next subsection. In 
relation to the insulating state, it is also an interesting 
future problem how the 120 degree Neel order, which is 
realized for U/t — > oo and t'/t 1, depends on U/t and 

6.2 as an RVB state 

Although \1/q has succeeded in describing a Mott tran- 
sition for an intermediate coupling strength, some adjust- 
ments in the spin correlation seem to be needed for the 
insulating regime {U > Uc), where a strong-correlation 
scheme is convenient. 

As a previous study for the square lattice showed, 
the strong-coupling regime {U > Uc for half filling) of 
the Hubbard model shares the properties with t-J-type 
models; the physics adiabatically connected to the limit 
of U/t oo. This picture also holds in the insulating 
state here, because the energies are scaled with t'^/U, 
and various quantities obtained with ^'g changes very 
slowly as a function of U/t and t'/t, as observed in n(k) 
[Fig. 5], N{q) [Fig. 7], D [Fig. 8] and ^(q) [Fig. 10]. In 
this sense, is substantially equivalent to Anderson's 
RVB state.ss 

For U/t oo, the present Hubbard model eq. (1) is 
reduced to the J- J' Heisenberg model with the same con- 
nectivity. Some approximations®^ to this spin model con- 
cluded that the Neel order persists up to J' / J ~ 0.6-0.7 
[t'/t = 0.77-0.84) beyond the classical value J'/J = 0.5 
{t'/t = 0.71). In contrast, in this study, the phase bound- 
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ary between the AF state and the insulating state of 
for ?7 ^ oo is situated at t' /t < 0.2, as deduced from 
Fig. 3(a). Thus, the two approaches quantitatively give 
different estimates. A possible origin of this discrepancy 
is that, despite the marked AF short-range correlation, 
'^Q does not have a seed of AF longer-range correlations. 
Also, the fact that the rcnormalization of t' /t is not intro- 
duced into ^fg^ possibly affects the discrepancy. As an 
improvement in both and ^q^, it is very important 
to study a wave function that can possess both rf-wave 
pairing and AF order. This will expand the area of AF 
insulator to some extent (see Fig. 14). 

Nonetheless, note that has an extremely low en- 
ergy as an insulating state; it readily defeats the AF 
state even with weak frustration, as in Fig. 3. In this 
connection, we need to remember Liang et a/.'s study, ^® 
in which the Heisenberg model on the square lattice is 
studied with intelligible wave functions, which are com- 
posed of a product of singlets with a bond-length dis- 
tribution as a variational parameter. They showed that 
the wave functions including long singlet bonds exhibit 
an AF long-range order with very low energies; the wave 
functions composed only of short-bond singlets do not 
have AF long-range order but have very close energies 
to the former. Thus, the stability of the wave function is 
little influenced by whether the spin correlation is short- 
ranged or long-ranged. There is no reason to suppose 
that this situation does not hold for frustrated cases. Be- 
cause Liang et al.^s wave functions are not tractable for 
t' /t 7^ 0, we should take an alternative approach, for ex- 
ample, refine a class of projected singlet wave functions 
by adjusting the gap parameters. They are expressed as 
the same form as eq. (7) with (/3k) which generally differs 
from that of the BCS type, eq. (8), but is tractable in 
the present scheme. 

7. Summary 

We have studied a Hubbard model on an anisotropic 
triangular lattice at half filling by an optimization VMC 
method with adequate precision, taking account of the 
SC and Mott transition arising in k-ET salts. In the 
trial functions — normal, AF, and d-wave singlet states, 
we introduc;e intersite correlation factors that control the 
binding between a doublon and a holon. We have suc- 
ceeded in describing the c?-wave SC and a Mott transition 
simultaneously in a sole approach, and in accounting for 
the behavior of a series of k-ET salts broadly. We itemize 
our main results: 

(1) Within the rf-wave singlet state, a first-order Mott 
(conductor-to-nonmagnetic insulator) transition takes 
place at Uc approximately of the bandwidth for arbitrary 
values of t' /t. The critical value Uc/t gradually increases 
as band width (or t' /t) increases. This transition is not 
directly related to a magnetic order. 

(2) The nonmagnetic insulating state (rf-wavc singlet 
state for U > Uc) has a considerably low energy and a 
strong short-range AF correlation. This state has a spin 
gap, as well as a charge gap. 

(3) The AF long-range order prevailing in the weakly 
frustrated cases {t' /t ^ 0.4) is rapidly destabilized as t' /t 
increases, and yields to the d-wave singlet state (SC or 



nonmagnetic insulating). 

(4) SC with the dx2_y2-wa.ve symmetry appears for 
moderate values of U/t (~ 6) and t'/t (0.4 ^ t'/t ^ 0.7), 
whose area is adjacent to both areas of an AF long-range 
order and a Mott insulator. The phase diagram obtained 
in this study is shown in Fig. 14. 

(5) SC pairs are formed with the aid of the short-range 
AF spin correlation, which is weakened by the frustration 
and vanishes for t'/t ^ 0.8. The SC transition is induced 
by the gain in interaction energy; this mechanism is iden- 
tical to that in the weak correlation limit as well as that 
of conventional BCS superconductors. 

In this paper, we have focused on the dx2_y2-wa.ve 
symmetry, which is established for t'/t = 0.*"° It has 
been reported, however, that the symmetry of pairing 
switches to d+{i)d-type symmetries at t'/t = 1.23,71,72 
It is hence worthwhile to check other pairing symmetries 
for large values of t'/t. Moreover, it is interesting to check 
with VMC how the features revealed in the single-band 
Hubbard model are altered in more realistic multiband 
models, especially in strong and intermediate correla- 
tion regimes. In connection with pairing symmetry, we 
finally comment that it is interesting to study the phase 
sensitive property, e.g., the tunneling effect and Joseph- 
son effect, in unconventional superconductor near a Mott 
transition.-'^'*' ''^ 

Acknowledgments 

The authors thank M. Ogata for useful discussions. 
They appreciate the discussions with and comments 
of H. Fukuyama, N. Nagaosa, K. Kanoda, K. Kuroki, 
and H. Kontani. This work is partly supported by 
Grant-in-Aids from the Ministry of Education, Cul- 
ture, Sports, Science and Technology, by NAREGI 
Nanoscience Project of the same Ministry, by the Su- 
percomputer Center, ISSP, University of Tokyo, and by 
a Grant-in- Aid for the 21st Century COE "Frontiers of 
Computational Science" . 

Appendix A: Gutzwiller approximation for d- 
wave state 

In this Appendix, we compare the results of conductor- 
insulator transitions in recent studies'^"' '""^ using GA for 
the d-wave BCS states^^ with the previous ones using GA 
for the Fermi sea. Two groups have recently addressed 
the present subject by applying GA for d-wavc BCS 
states to Hubbard-type models with extra exchange- 
coupling terms. These extra terms are added for tech- 
nical reasons, and seem to make no essential difference. 
Powell and McKenzie^* studied a Mott transition and 
SC without allowing for AF. They found a first-order 
SC-insidator transition at U/t somewhat larger than the 
bandwidth. In the insulating regime, however, the kinetic 
energy vanishes, indicating this transition belongs to the 
Brinkman-Rice type. The AF effect is incorporated into 
GA for the d-wave BCS state by Gan et al.,^° who found 
a SC-AF insulator transition at U/t somewhat smaller 
than the bandwidth. In the SC regime, the SC correla- 
tion and rcnormalization of g develop as U/t increases 
without sublattice magnetization, whereas in the AF in- 
sulating regime, there appears a finite sublattice magne- 
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tization but g is reduced to the value at C//t = 0, namely 
g = \.ln fact, this result is closely parallel to that of the 
GA for the Fermi sea into which an AF effect is intro- 
duced, as in the following. 

About three decades ago, Ogawa et alJ^ introduced 
the AF effect into the GA formula for the Fermi sea, us- 
ing ■pQ^AF, with a cubic-type lattice in mind. The result 
obtained exhibits a transition; for < C/ < C/afj the re- 
sult coincides with the ordinary (or paramagnetic) GA 
formula (Aaf = 0, g < 1), whereas for U > Uaf, the 
result is switched to the mean-field solution of AF insu- 
lator (Aaf > 0, g = 1). Here, Uaf denotes U at which 
the energy curve of the GA formula intersects with that 
of the AF mean-field approximation. Thus, the transition 
at ?7 = Uaf is a first-order AF transition which does not 
correspond to the Mott transition we considered in this 
paper, but is basically identical to those derived by the 
mean-field theories for t' ^ 0.^^ Later, it is shown by a 
VMC study''^ for Vq^af that this conclusion is again 
incorrect; the true variational state has an AF gap and 
a finite renormalization by g (Aaf > 0, 5 < 1) simul- 
taneously for arbitrary values of U{> 0), and does not 
exhibit a transition for U > 0. 

Thus, as long as a Mott transition is concerned, the 
results of GA for the d-wave BCS state''"- share the 
features with the conventional GA for the Fermi sea.""'' ''^ 
We should be sufficiently careful to use GA to study a 
Mott transition. Nevertheless, it is fair to mention that 
the results of GA with respect to SC are often consistent 
with the present ones shown in §5. 

Appendix B: Mott transition by doublon-holon 
binding factors 

Up to this time, a Mott transition has been repeatedly 
studied using the doublon-holon binding factor Vq and 
its minor variations. However, the course of study has 
not been straight but a little confusing. In this Appendix, 
we first recapitulate this subject (§B.l), and then make 
comments on misunderstandings based on the above con- 
fusion in a recent VMC study^^ (§B.2). 

B. 1 Mott transition by Vq 

For simplicity, let us consider the case oi t' — [ji' — 
0), which makes no substantial difference for finite t' . 
Because in eq. (4) is close to 1 near a Mott transition, 
the wave function eq. (6) is expanded with respect 
to 1 — /z (= m) as 



(0) _ 
Q - 



(1) _ 
Q - 



[n(i-o^)]*G, 

i 



(B-1) 
(B-2) 

(B-3) 



(2) 

Q 



= { E n (i-Qi)]}*G. (B-4) 



The bases that compose possess necessarily £ lattice 
sites which are occupied by doublons (or holons) without 
a holon (doublon) in their nearest-neighbor sites. Note 



that in the vicinity of a Mott transition, g (Gutzwiller 
factor) is so small that the doublon density is very low 
{D ^ 0.04, as shown in Fig. 8). Using eq. (B-1), the 
kinetic and interaction energies are formally expanded 
with respect to to, as 

Kq + mKi + TO^ifa + • • • , (B-5) 

/o + m^h + ■■■ , 



E\iin{g,m) 
Eint{g,m) 



(B-6) 



in which Kg and {£ = 0,1,2, ■■ ■) are functions of the 
Gutzwiller factor g 

Ko + Io=Eo = {0\n\0)/No, (B-7) 

Ki = [(0|Hki„|l) + (l|Hki„|0)] /No, (B-8) 

K2+l2= [(0|Wkin|2) + (l|Wkin|l) + (2|Wkin|0) + 



(l|Hi„t|l)+iVi(0|H|0) /No 



(B-9) 



with |*g^) = \i) and Ne = Thus, the total energy 
is given as 

E{g, m)^Eo + mKi + {K2 + h) + ■ ■ ■ ■ (B-10) 

From eq. (B-10), it seems that E{g,m) has a mini- 
mum at a finite m (or < 1), if Jfi is negative. Ac- 
tually, Yokoyama and Shiba^^ calculated Ki for one- 
dimensional systems, and confirmed that Ki is negative. 
They thus concluded that m is always finite as long as 
U /t is finite, that is, a Mott transition does not arise. 

Later, however, this conclusion turned out to be sub- 
stantially incorrect through a more careful study. Be- 
cause, at half filling, a doublon and a holon arc cre- 
ated simultaneously as a pair, a basis included in 



(1) 
Q 



( 2 ) 

has at least two doublons, whereas \E'q contains bases 
with only one doublon, as shown in Fig. B-1. Thus, as 
g decreases, \Ki\ can become negligible, compared with 
\K2\- In fact, Ki/t has been calculated using VMC with 
VqVg^f for the one-dimensional Hubbard model, with 
the result that as g decreases, \Ki\/t is abruptly reduced 
from ~ 0.1 (5 = 0.4) to ~ 0.007 {g = 0.2).2i In addi- 
tion, the quadratic behavior of E{m,g) for small g and 
m ^ 0, namely dE/dm = 0, is clearly observed for two- 
dimensional cases. Eventually, it was confirmed that 
Mott transitions actually take place for the square lat- 
tice Hubbard model by applying Vq to ^p^" and ^d-^'^ 
Afterward, Mott transitions have been observed using 
Vq for various systems in addition to the present model: 
the lattice shown in Fig. 1(b) with both $f and ^d,^^ the 
kagome lattice with $f,'^* the checker-board lattice with 
$F,^^ and a degenerate Hubbard model on the square 
lattice with ^f-**" 

Independently of ref. 21, Millis and Coppersmith^^ 
studied the possibility of a Mott transition with a vir- 
tually equivalent variational wave function to VqVg'I^f- 
To determine whether or not a Mott transition arises, 
they adopted a general theory of Kohn^^ in estimating a 
zero-frequency part of the optical conductivity from the 
linear response of twisted-boundary systems to a small 
flux. They concluded that a wave function like VqVq^f 
never becomes insulating for a finite U/t, on the basis of 
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Fig. B-1. Examples of electronic configurations appearing (a) in 



and (b) in "^q' , at half filling. Solid (open) circles indicate 
doublons (holons). Each example belongs to the case in which 
the doublon number is minimum. 
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Fig. B-2. (Color online) Comparison of VMC results for t'/t = 
0.7 between the study of Liu et al. (taken from Fig. 1 in ref. 23) 
and the present study for (a) the d-wave gap (the optimized 
variational parameters), and (b) the pairing correlation function. 
The quantities studied by Liu et al., Ao/t* and ips, (crosses) 
practically correspond to A^/t and PJ^™ in the present study 
(solid symbols), respectively. 



a VMC result with a nonoptimized g{— 0.3) for a one- 
dimensional small (10-site) system. However, this con- 
clusion is an error, as we have demonstrated with actual 
counterevidence in the preceding paragraph (or in §4). 
We need to check this study as to which process causes 
the error. 

B.2 Comments on recent VMC study 

Recently, Liu et al.^^ have studied the same model 
eq. (1) using a VMC method with an almost identical d- 
wave singlet trial function (abbreviated as ^'liu) to ours, 
\1/q, in eq. (6). The differences between \1/Liu and are 
as follows: (i) In ^'liu, to improve the renormalization 
effect of the quasi-Fermi surface, t" is additionally intro- 
duced in the other diagonal (—1,1) direction in Ek as a 
parameter, (ii) In ^'liu, the magnitude of the doublon- 
holon binding factors for next-nearest-neighbor sites is 
assumed identical to that for the nearest-neighbor ones. 



namely /i' = /i in our notation. Because these differ- 
ences are not essential,™ it is natural that the results of 
Liu et al. are not profoundly different from the present 
ones, and that some of their arguments are common to 
ours; we are ready to respect this aspect of their study. 
However, in a crucial point in connection with a Mott 
transition, their interpretation is distinct from ours and 
questionable. Because it is vital to clarify the properties 
of Mott transitions in this type of wave function for re- 
lated studies, it is necessary to clarify the questionable 
points. 

In the study of Liu et al., a pair correlation function 
(ips) was calculated for t'/t = 0.7 with a fixed system 
size, L = 12. We have traced 4's from Fig. 1 in ref. 23, 
and plotted it in Fig. B-2(b) together with our corre- 
sponding results. For small values of U /t, tps is negligible, 
whereas it increases at approximately U/t = 8.5 and then 
it gradually decreases. This behavior of ips is interpreted 
in ref. 23 as follows: For small values of U/t, the state is 
a conventional metal; robust SC occurs at approximately 
U/t = 8.5, and weak SC continues to larger values of U/t 
without a transition. The problematic interpretation con- 
sists in the last point. In comparing tps with P^™, U/t at 
which Ips becomes conspicuously large (~ 8.5) is larger 
roughly by 1 than Uc/t (=7.45) of PJ'™. This discrep- 
ancy may stem from point (ii) of the differences in the 
wave functions, and is of little importance. The most 
crucial point is that they overlooked the Mott transition 
this type of wave functions exhibits immediately above 
the region of robust SC. They probably accepted the in- 
correct claim of Millis and Coppersmith mentioned in 
§B.l, and neglected to check the system-size dependence 
of ■(/'s.^° For such a finite system as L = 12, tps is still 
visibly finite, but must decrease abruptly with increas- 
ing L, as we showed the system-size dependence of PJ'™ 
in the inset of Fig. B-2(b). 

This misinterpretation influences some other points. 
Because the region of SC is considered to continue to 
larger values of U/t ( ^ 8.7), the features of the insu- 
lating state are confused with those of the SC state. Al- 
though the prominent renormalization of Sk (or f and t") 
is considered to occurs in the SC state, as shown in Fig. 3 
of ref. 23; correctly, the renormalization of t'/t is small 
in the SC state, but prominent in the insulating state, 
as we showed in §4 and §5. Furthermore, the spin corre- 
lation function is compared between the metallic regime 
{U/t = 7), for which tps is negligibly small, and the insu- 
lating regime {U/t — 11) in Fig. 4 of ref. 23. Therefore, as 
long as the spin correlation in the SC state is concerned, 
this comparison is useless. 

Finally, we point out the possible fictitious enhance- 
ment of PJ'™ (or Ips) specific to finite system sizes. In 
Fig. B-2(b), the values of P^™ immediately below Uc for 
L = 10 are roughly twice larger than those for L = 12 
and 14. This is because the A;-point configuration for 
L = 10 happens to be specially advantageous for SC un- 
der the given condition, including the renormalized effect 
of ek by U. This is confirmed by the fact that the quasi- 
Fermi surface has a different shape only in these special 
cases (not shown). Undoubtedly, the behavior of PJ^™ for 
L = 12 and 14 will be closer to that for L = oo. In accor- 
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dance with this anomalous enhancement in P^™, A^/t 
takes large values as observed in the inset of Fig. B-2(a). 
According to our experience, such behavior is not rare 
for the marginal values of t' jt where SC is vanishing. 
When we return to the data of in Fig- B-2(a) with 
this situation in mind, we notice that the values of Aq/I* 
for U jt ^ 8.5 (three points), where; (.i's is considerably en- 
hanced, are specifically increased. Hence, the prominent 
peak in ■0s at iJ/t ~ 8.5 is probably a fictitious increase 
specific to their system size and boundary conditions. 
We hope improved studies of will appear. 
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